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1. Introduction 

High energy polarized e + e~ colliders will be essential instruments in the search for the fun- 
damental constituents of matter and their interactions. One such collider being designed 
is the International Linear Collider (ILC) which is expected to run at centre-of-mass en- 
ergies > 500 GeV. At these energies there will be a significant proportion of top quark 
pairs produced from the annihilation process so the ILC provides an impressive tool to 
carry out detailed studies of top quark physics. The top quark radiates gluons both in its 
production phase and its decay phase, thus it is useful to see how the leading order ex- 
perimental analysis will be affected by these QCD corrections. There have been numerous 
studies of top quark production and their decays both at leading order and next-to leading 
order in QCD. Details can be found in [1-8] and many more besides. In this paper we 
consider the process at next-to-leading order in the production and semi-leptonic decays 
of the top pairs using the POWHEG method [9, 10] in conjunction with the Monte Carlo 
event generator Herwig++ [11]. The POWHEG method has been successfully applied to Z 
pair production [12], heavy flavour production [13], e + e~ annihilation into hadrons and 
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Drell-Yan vector boson production [14,15]. We work in the narrow width approximation 
and hence do not include interference between the production and decay emissions which 
are negligible in this limit [16, 17]. We also take account of the beam polarization and spin 
correlations of the top pairs. Finally, we present plots of some relevant distributions. 



2. Hardest emission generation: Production 



The order-a s differential cross section for the process e + e 
a vector current, can be written as 



V — ► ttg where V represents 



R(x,y) = ayWy(x,y) 



ay 2a s 
v 3ir 



(x + 2p) 2 + (y + 2pf + C v 2p 
(l + 2p)(l-z)(l-y) (1 



2p 



x) (1-yY 

(2.1) 

where ay is the Born cross section for heavy quark production, Wy = R/ay, p = m t 2 /s 
where m t is the mass of the top quark and s is the square of the center of mass energy, 
(y = — 8p(l + 2p), v = y/1 — 4p and x, y are the energy fractions of t and i respectively. 
In the case of the axial current contribution e + e _ — > A — ► tig, we have 



R(x,y) = a A W A {x,y) 



a A 2a s 
v 3ir 



(x + 2pf + (y + 2pf + ( A 



2p 



2p 



(1-4^(1-^(1-2/) (l- X f (l-yf 



(2.2) 

where a a is the Born cross section for heavy quark production by the axial current, Wa = 
R/aA and (a = 2p[(3 + x g ) 2 — 19 + 4p] where x g = 2 — x — y. 

Because of the top mass, the phase space for gluon emission is reduced and the collinear 
divergences present in the massless quark cross-section are regularized here. However, the 
infra-red divergence as the gluon momentum goes to zero is still present. We can write 
down the cross-section for the hardest emission as 



da = J2 B ( v ) d ®v 



Af LO \0) + A% 



(NLO) 



R(v,r) 



(2.3) 



where B(v) is the Born cross section and v represents the Born variables, r represents 
the radiation variables and d$ v and d& r are the Born and real emission phase spaces 
respectively. A^ lo (pt) is the modified Sudakov form factor for the hardest emission with 
transverse momentum pr, as indicated by the Heaviside function in the exponent of (2.4), 



a NLO 



(pr) = exp 



d$r^^-@(k T (v,r)-p T ) 
B{v) 



Furthermore, 



B{v) = B{v) + V(v) + / (R(v, r) - C(v, r))d$ 



(2.4) 



(2.5) 



B{v) is the sum of the Born, B(v), virtual, V{v) and real, R(v,r) terms, (with some 
counter-terms, C(v,r)). It overcomes the problem of negative weights since in the region 
where B(v) is negative, the NLO negative terms must have overcome the Born term and 
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hence perturbation theory must have failed. It is used to generate the variables of the Born 
subprocess to which the real-emission contributions factorize in the collinear limit. 
Now explicitly for e + e~ — > tig, 



K L °(pt) = exp 



dx dyW(x, y)@(k T (x, y)) - p T ) 



where 



kr(x,y) = \ s 



(l-x)(l-y)(x + y-l)-p{2-x- yf 



max(x, y) 2 — Ap 



(2.6) 



(2.7) 



is the transverse momentum of the hardest emitted gluon relative to the splitting axis, as 
illustrated in Figure 1 below. 




Figure 1: Transverse momentum, Ut- 



2.1 Generation of radiation variables, x and y 

The radiation variables, x and y are to be generated according to the probability distribu- 
tion 

A w (k T )W(x,y)dxdy (2.8) 

where for e + e~ annihilation via a vector and axial current, W(x, y) and A w '(kr) are defined 
in (2.1), (2.2) and (2.6). 

In the region where x > y, let us define the dimensionless variable, k as 



K 



k 2 T (1 - x)(l - y)(x + y - 1) - p(2 - x - yf 



S X 2 — 4/9 

There are two solutions for y for each value of x and k. 



(2.9) 



2/1,2 



x 2 - 3x - 2px + 2 + 4p ± yj{x 2 - 4/))(4k(x - 1 - p) + (x - l) 2 ) 

2(1 + p-x) 



(2.10) 
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Exchanging the y variable for k, we find 

J W(x,y)e(k T (x,y)- PT )dxdy = J*™ dx £™ dn^^-^W(x, y h2 ) 

/f , 2a s (Ks) \/x 2 — 4p . . 

dx / dK iL j V / ==W x,y 1>2 

J 37T y^/^X - 1 - p) - (x - l) 2 ) 

(2.11) 

Now the integrand, W' = V x ~ 4 p = W(x,yi2) in (2-11) yields a complicated 

integral so we look for an upper bound on W' which we denote as V' = ^V(x, 2/1,2) to 
simplify the integration. We then generate the radiation variables as outlined below: 

1. Set p max = fcrmax- 



2. For a random number, n between and 1, solve the equation below for px 

where A v (pr) = exp — j dxdnV ' (x, k) 

3. Generate the variables x and y according to the distribution 

V(x,y)5(k T (x,y)- PT ) . (2.13) 

4. Accept the generated value of pt with probability W/V. If the event is rejected set 
Pmax = Pt and go to 2). 

Using our knowledge of the form of the integrand in the massless quark case [18], we 
guess that V' should take the form, 

T r 1 , \ « T 2a s (ns) 4 .„ _ 

V ») = JV. d l (1 _ i + 7(K i))7()tii . ) (2.14) 

where 

7(k, x) = ^/(l - x)(l - 2/c - 2v^2 + /9/e) . (2.15) 

and is a normalisation factor which depends on k, which has to be tuned to ensure 
that V' is an upper bound of W' . Both V' and W' have the same divergent behaviour at 

Imax = 1 — 2/C — 2y / K 2 +~pK. 

Table 1 shows the 7V K values for different ranges of k and the two solutions for y used for 
both axial and vector currents. The lower limit on n was set by choosing kr = Aqcd = 0.2 
GeV thus setting a lower bound on the transverse momentum. 

We now consider the specific case where nit = 175 GeV and -^/s = 500 GeV i.e. 
p = 0.1225. For k < 0.024, there are two y solutions in the region of phase space where 
x > y. This is illustrated in Figure 2 below for k = 0.01. The red line denotes the phase 
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Range of k 


AT (nu \ 


AT fn„\ 


0.024 - 0.03 


0.4 


0.4 


0.015 - 0.024 


0.7 


0.7 


n 005 — 01 5 


1 9 


i i 

_1_ . _1_ 


0.0005 - 0.005 


4.0 


2.6 


0.0001 - 0.0005 


9.0 


6.0 


0.00005 - 0.0001 


13.0 


7.0 


0.00003 - 0.00005 


17.0 


10.0 


0.00000016 - 0.00003 


45.0 


35.0 



Table 1: N K for different values of k for both axial and vector currents 




Figure 2: Phase space and y solutions for k = 0.01 in the region x > y. 



space for gluon emissions. The two solutions lie on either side of the dashed line 

y= (g-;0(i-' + 2p) (2 . 16) 

y 2(1- x + p) K J 

and are equal when 



x 



1 -2k- 2yjn 2 + pK (2.17) 



which lies on the dotted line. At k = 0.024, the branches meet along the line y = x and 
there is only one solution for y in the region (the lower branch). So for k > 0.024, only 
one y solution exists for x > y. This is illustrated in Figure 3 for k = 0.028 . In addition 
there are no y solutions for n > 0.03 for x > y . Also note that for x < x u , there is only 
one solution for y. 
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x 



Figure 3: Phase space and y solutions for k = 0.028 in the region x > y. 

In the region where there are two solutions, the integral (with V' in place of W ) in 
(2.11) is performed along both branches independently and summed. For the upper branch, 
x runs from x u to x max = 1 — 2k — 2a/ k 2 + pn while for the lower branch, x runs from x l 
to Xmax where if we define 



x a = 39k - 1 + 12/9 - 168/5K - 48p 2 + k 3 + 15k 2 + 12 P k 2 + A8p 2 K + 64p 3 , 
x b = 6(-33k 2 + 3k - 3k 3 + 288p 2 K - 768p 3 K - 48/w + 168pK 3 - 204p 2 K 2 + 300pK 
+ 12/9K 4 + 768 p 4 K + 1U P 2 k 3 + 486p 3 K 2 )^ , 



x c = V 'x a 2 + x 6 

, -1 f x b 
Xd = tan I — 



x f 



1 \ (X d ^ 
Xr COS 

12 V 3 

(-1 -k 2 - 10k + 8p - 2 P k - IQp 2 ) cos (^) 

T ' 
12x c 3 

re + 5 + Ap 
6 ' 

= t - m (xi + 1+^+1^-^+2^+16^ (218) 



12 V 3 
we can write x u and as 



X Xg "T" Xy ~\~ Xg ~t~ X/j_ 5 

x l = x e + Xf + Xg — Xh] (2-19) 
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In the region where there is only one solution for y, x runs from x l to x u . 
The k integration can then be performed numerically. Having performed the integration, 
values for k and hence kx are then generated according to steps 1 and 2 in Section 2.1. 
The variables x and y are then to be distributed according to W'(x, y)S(kT(x, y) — Pt)- 
This is the subject of the next section. 

2.2 Distributing x and y according to W(x, y) 

To generate x and y values with a distribution proportional to V(x,y)8(kT{x,y) — Pt), 
where from (2.11), V{x,y) is the d yj dn , we can use the S function to eliminate the y 
variable by computing 

V(x,y) 



D(x) = J dy5(k T -p T )V(x,y) = 



dk T 



(2.20) 



y=y 

where y is such that kr{x,y) = pr- Note that is the same for both y solutions. We 
then generate x values with a probability distribution proportional to D with hit-and-miss 
techniques as described below. All events generated have uniform weights. 

1. Randomly sample x, N x times (we used N x = 10 5 ) in the range [x m i n : x max ] for the 
selected value of k. 



2. For each value of x, evaluate D = D(x,y±) + D(x, y-i) if there are 2 solutions for the 
selected k and D = D(x,y2) if there is only one solution. Also, if k < 0.024 and 
x < x u (see Figure 2), there is only one y solution so evaluate D = D(x,y2)- 

3. Find the maximum value -D max of D for the selected value of n. 

4. Next, select a value for x in the allowed range and evaluate D. 

5. If D > rD max (where r is a random number between and 1), accept the event, 
otherwise go to 4.) and generate a new value for x. 

6. If for the chosen value of x, there are two solutions for y, select a value for y in the 
ratio D(x,yi) : D(x,y 2 ). 

7. Compare V' (x,y) with the true integrand, W'(x,y). If the event fails this veto, set 
n max = k and regenerate a new k value as discussed in Section 2.1. 



NB: For the region y > x, exchange x and y in the above discussion. In this way, the 
smooth phase space distribution in Figure 5 below was obtained for the hardest emission 
events for an axial current. The plot show 2,500 of these events. 
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Figure 4: Phase space and distribution of hardest emissions for axial (left) and vector(right) 
currents with p = 0.1225. 




+ + * ++ + t + + +++ + +t # # 




Figure 5: Phase space and distribution of hardest emissions for an axial current with p 
(left) and p = 0.01361 (right). 



0.0625 



The procedure was repeated for y/s = 700 GeV, p = 0.0625 and ^fs = 1500 GeV, 
p = 0.01361 and the corresponding plots are shown below. As can be seen, the method is 
stable as p — ► 0. This is not surprising because the upper bound function V(x, k) in (2.14) 
is stable as p — > and tends to W' , the true value of the integrand in this limit. 
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3. Hardest emission generation: Decays 



In addition to gluon emission in top production, we also studied the emission in its decay, 

t(pi) - W + ( Wl )b( ri )g(k) . (3.1) 

The procedure for generating the hardest emission in this case follows the same lines as 
discussed in Section 2. We parameterize the phase space for the decay in terms of variables 
x and y defined as 

2wi • p\ 

y = — 2 — a 

mf 

x = ^ (3.2) 

m t 

where a = m 2 w jm\ with m w and mt the masses of the W boson and top quark respectively. 
(y + a)/2 and x/2 are the energy fractions of the W boson and gluon in the top frame. 
Therefore the corresponding energy fraction of the b quark in this frame is given by 

Xh 2 — y — a — x 

i = — 2 — • (3 - 3) 

In this paper, we neglect the b mass and work in the narrow-width approximation so that 
the top quarks and W boson are on-shell. The t — ► Wbg differential decay rate is given by: 



1 d 2 T a s C F 



(1 - y)(l - x) +x 2 (y + x-1) 2 2a(l-y)x 2 



T dxdy ir(l-y)x 2 [ I- a 2(1 - a) 2 (l-a) 2 (l + 2a) 

(3.4) 

where Tq is the leading order decay rate. The phase space limits for the decay are: 

(XX , » 

+ (1 - x) < y < 1 , 



1 — x 

< x < 1 - a . (3.5) 

Working in the rest frame of the top quark where the parton shower is formulated in 
Herwig++, we identify the splitting axis corresponding to the original b — W boson axis 
and therefore the relative transverse momentum for gluon emission is: 



fc / \ 0- - V)(y + x(2 - y - a) - x 2 - 1) 

hr(x,y) =m t y {y + a ? - 4a " (3 ' 6) 

k 2 

Now defining a dimensionless variable k = we find that in analogy to the production 
case, there are 2 solutions for y for each value of x and k. 



_ x 2 + ax + 2 - 3x - 2qk ± ^(x 2 - 4k(1 + a))(x - l) 2 + 4c«;(4k + 1 - a) + x 2 (a + 2x - 2) 
2/1,2 ~ 2{k + 1-x) 

(3.7) 

These solutions may be identified with either initial state gluon emission from the top 
quark (3/2) or final state radiation from the bottom quark (yi). A plot of the phase space 
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Figure 6: Phase space(solid), y (dot-dash) and solutions y\ (dots) and y2 (dashes). 



and the 2 solutions for k = 0.01 is shown in Figure 6. We then construct the modified 
Sudakov form factor for the generation of the hardest emission. The exponent of the form 
factor is given by 

f W(x,y)G(k T (x,y)- PT )dxdy = dx d ^^^ ^W(x, k) , (3.8) 

J Jx min Jk IT dK 

where W(x, k) is the differential cross-section (3.4) and ^ is the Jacobian for the change 
of variables from y to k. Note that k = K max = — when the W boson is at rest and 
x = 1 — ^/a, y = 2yfa. For a given k, x m \ n and x max are also given by 



^min 
^max 



2v^ 

1 — a — 2\J7m . 



(3.9) 



To make the integral simpler, we again look for an upper bound V (x, k) on the integrand 
as we did for the production case. To do this we replace the Jacobian with the simpler 
expression, 

dy d ( x 2 — 3x — 2na + 2 + xa\ —a x 2 — 3x — 2na + 2 + xa 

~d^ = d^\ 2(k + 1-x) ) = k + 1-x 2(k + 1-x) 2 ' ( ' ' 

where y lies in between y\ and y<i and is indicated in Figure 6. We also overestimate the 
differential cross-section by replacing (3.4) with 

olsCf 1 — a 



U(x,y)=N K - 



it 2x 2 (l — y 



(3.11) 



where N K is a normalisation factor dependent on k and is chosen such that V' = U^- is 
greater than the integrand in (3.8). The N K values are given in Table 2 for the 2 solutions. 
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Range of k 



N K { yi ) x 10 4 N K {y 2 ) x 10 4 



0.01 - 0.0737 


0.005 


0.006 


0.005 - 0.01 


0.0175 


0.02 


0.001 - 0.005 


0.03 


0.045 


0.0001 - 0.001 


0.08 


0.12 


0.00005 - 0.0001 


0.2 


0.2 


0.000025 - 0.00005 


0.3 


0.2 


0.0000075 - 0.000025 


1.0 


0.9 


0.000005 - 0.0000075 


2.0 


0.9 


0.0000025 - 0.000005 


3.0 


0.9 


0.0000013 - 0.0000025 


6.0 


0.9 



Table 2: N K for different values of n 
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Figure 7: Phase space distribution of POWHEG events 



The lower limit on k is 1.3 x 10~ 6 and was set by choosing kx = Aqcd = 0.2 GeV thus 
setting a lower bound on the transverse momentum. We then generate the values of k and 
distribute x and y according to the true differential (3.4) using vetoes as described for the 
production case in Sections 2. Figure 7 shows the phase space distribution obtained. 

4. Spin Correlations and the distribution of Born variables 

In [19], it was observed that the lepton matrix element for the production process 

e + (p) +e~{q) -► t( Pl ) + i(p 2 ) + g(p 3 ) -»• W + ( Wl ) + 6(n) + W-{w 2 ) + b{r 2 ) + g(p 3 ) 

-»• l + (h) + v(x x ) + 6(n) + r(k 2 ) + u{x 2 ) + b{r 2 ) + g(p 3 ) (4.1) 
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is bounded from above in the narrow width approximation by the undecayed matrix ele- 
ment obtained by eliminating the decay products i.e. W + , W~ , b, b and putting the parent 
particles i.e. t,i on-shell, multiplied by a process dependent constant. We can then use the 
undecayed matrix elements to perform computer-intensive tasks such as event generation 
and finally, by using the hit-and-miss method, replace the parent particles with their decay 
products. This procedure is outlined below: 

1. Evaluate the undecayed matrix elements which are proportional to the upper bounds 
on the lepton matrix elements. Generate hard events using the POWHEG method 
described above with the top and anti-top quarks in the final state. 

2. For each event, generate the decay products and their four-momenta according to 
the phase space. 

3. Evaluate the leptonic decay matrix element for each event. If the decay matrix 
element divided by the corresponding upper bound is less than a random number r 
between and 1, throw away the decay momenta and return to step 2. 

4. Otherwise, replace the top and anti-top momenta with the decay momenta and shower 
the event. 



4.1 Undecayed matrix elements 

At the ILC, the electron and positron beams will be polarized i.e. either e^e^ or e^ej 
where the subscripts L and R represent the left-handed and right-handed helicity states 
respectively. The corresponding undecayed matrix element for e£e^ annihilation is: 



[v{q)l^Lu{p)]u{pi,s t ) 



I — —1 L + -TT^R J (-P2 -P3 + mth„ 



2p 2 • p 3 v 2 



(O-LL n dLR 



+ 7^ 7v(Pi + P3 + mt) -ttIl + 

2p 1 -p 3 v 2 



T a v{p 2 ,st)el{p 3 ) 



(4.2) 



where e is the polarization vector of the gluon, T a is the colour matrix, p = p^j^ and 
7r/l = 7^(1 ^ 75)/2- For e^e^ annihilation, interchange L,R in the above equation. st,Sf 
are the spin vectors of the top and anti-top quarks respectively and satisfy the relations: 



st - Pi 

St ■ P2 

St ■ s t 
St ■ St 



(4.3) 
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The massive spinors u(p, s),v(p, s) are given in terms of the massless spinors u(p),v(p) by 



u(p, T) 

u(p, |) 

v{pA) 

v(p,[) 



2 «(P) 

1-75S , . 
2 «(P) 

1 + 75S , , 
2 V (P) 

1 -75g , , 
2 V (P) 



(4.4) 



The coupling constants ajj are given by 



e 2 g 

au = 

s 



-Qt + QlQt- 



sin 2 6 W s-M 2 z + iM z T z _ 



(4.5) 



where Mz is the Z boson mass, Tz is the width of the Z boson, 9w is the Weinberg angle, 
Qt is the electric charge of the top in units of the electric charge e, g = y/Anas and s is 
the center of mass energy squared. The couplings to the Z boson are given by 



Qe 

Qt 
Qf 



2 sin 2 9 W - 1 
2 cos 6w 

sin 2 6\v 

cos Ow 

3-4sin 2 9 W 
6 cos 6w 

2 sin 2 9 W 

3 cos Ow 



(4.6) 



In Section 2.2, we distributed our events according to the vector and axial vector 
current matrix elements separately using the POWHEG method. To obtain a full unpolar- 
ized distribution we can select events from either current distribution according to their 
contributions to the full cross-section given below. 

a = 3/3(1 + 2p) + ci^) a VV + 3/3 3 (l + a AA 
47ra 2 



<?AA 



Anal 



[Ql - 2Q t V e V tX i(s) + (A 2 e + V e 2 )V t 2 X2 (s)] 
[(A 2 e + V 2 )A 2 tX2 (s)] , 



(4.7) 



where (3 = ^1 — p, a em is the electromagnetic coupling, At, A e and Vt, V e are the axial and 
vector coupling constants of the top t and electron e to the Z boson and c\ = 3.5 and 
d\ = 2.25 are the QCD correction coefficients defined at rat = 175 GeV and yfs = 500 GeV 
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i.e. p = 0.1225 [20]. Xi( s ) an d X2(s) are given by 



a(s - M 2 



ZJ 



(s - M\y + r|Mf 

X2(S) = "\s-M 2 ) 2 + T 2 z M 2 

«=^fi, (4.8) 
167ra e m 

where is the Fermi constant and and are the mass and decay width of the Z 
boson respectively. 

Explicit expressions for the Born, virtual and real polarization dependent squared 
matrix elements for the production process are given in [5] in terms of the energy fractions 
x, y of the top and anti-top quarks and the polar angle and azimuthal angles orienting 
the tig plane relative to the e + e~ beam axis. For each initial polarization, we then assign 
final-state polarizations to each event in proportion to the squared matrix elements and 
distribute the polar and azimuthal angles of the top/anti-top pairs accordingly using well- 
known Monte Carlo techniques. 

4.2 Decay matrix elements 

Next, we investigate the decays of the top and anti-top pair. The leptonic matrix elements 
for the process in (4.1) are dependent on the spins of the top and anti-top quark. This 
dependence can be written in the form of a decay density matrix. The decay density matrix 
p x y , for an on-shell top quark is given by 

4gM 



Jw y tb 

P\,X' ~ 



(wj - ml) 2 + (m w T w ) 2 

(ri • x 1 )(pi ■ ki) - (s t ■ fci)(ri • xi)m t -(k x ■ n)(x 1 ■ n)m t - ie(pi, fei, s t ,n)(x 1 ■ n) 
(hi ■ n){x\ ■ n)mt + ie(pi, ki,s t ,n)(xi ■ n) (n • xi)(pi ■ h) + (s t ■ fci)(n • xi)m t 

where A, A' are spin labels, st is the top spin vector, n is a spacelike vector perpendicular to 
st and pi and m w , T]y are the mass and width of the W boson respectively. In this paper 
we work in the helicity basis for which the top quark spin is defined along its direction 
of motion. A similar matrix can be derived for i decay. The spin-specific decay matrix 
elements are therefore of the form: 

5 W ;a;- = Mx t x,P XtX ' AlX 'Ml [K (4.9) 

where Aj,Af are spin labels for the top and anti-top respectively and M is the matrix 
element for the undecayed process introduced in Section 4.1. By diagonalizing the density 
matrix, we can obtain the largest possible value of the matrix elements and hence the upper 
bound. An explicit computation gives this upper bound | M* fe | 2 on the top decay as [19], 

I M \ \ 2 = 4^|T4 fe | 2 (ry-xiKgivfci) ^ 2 

1 ub 1 [(w 2 -m 2 w ) 2 + (m w T w ) 2 ] [(p 2 - m 2 ) 2 - (m t T t ) 2 ] 11 1 " j 
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where | M | 2 is the undecayed matrix element for unpolarized tig production. A similar 
expression | M l ub | 2 can be obtained for the decay of the top anti-quark by interchanging 
the labels 1 and 2 and t and t in (4.10). Hence the full upper bound can be written as: 

I M t l 2 l /If* I 2 

| M tt |2 = \M ub \\ M ub\ (4 n) 

| M | 2 

Having obtained the decay matrix elements and their upper bounds, we then proceed 
to generate events with leptons in the final state as outlined at the beginning of this section. 

For the PDWHEG decays, we apply the same method where in this case the undecayed 
matrix elements \M\ are the leading order matrix elements for the process 

e + + e~^t + i. (4.12) 



We then use the next-to-leading order decay matrix for which the helicity amplitudes can 
be found in [5]. These are given in terms of the polar and azimuthal angles of the decay 
w.r.t the top/anti-top axis and we distribute them as described for the production process 
in Section 4.1. Note that in this case, we generate two decay gluons, one each from the top 
and anti-top quark. 

In addition, we also consider POWHEG radiations in both the production and decay 
process by independently generating the emission and distributing the Born variables of 
the production process first and then generating the emission and distributing the Born 
variables of the decay process to yield three gluons in the final state. 



5. Decay NLO lepton spectra comparisons 



Extensive studies have been carried out on the lepton angular and energy distributions 
from the semi-leptonic decays of polarized top and anti-top quarks at next-to-leading order 
in a s [1,3]. 

t -> W + + b + g^e + + v e + b + g 

i W~ + b + g e" + v e + b + g . (5.1) 

In the top rest frame, we define 9 as the angle between the spin 3- vectors st,sj of the 
decaying quark and the lepton. We also defined the scaled energies x^ n of the charged 
lepton and the neutrino respectively as 

2E l 
xi = 

m t 

x n = 5.2) 

m t 

where E[ and E n are the energies of the charged lepton and neutrino in the top rest 
frame. In these variables, the NLO double differential distribution of the charged lepton 
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and neutrino in the decay of a heavy top or anti-top quark with polarization S has been 
shown to be of the form 



dr 



Ln 



dxi n d cos 6 



32vr 3 
2a s 



Fo' n (xi,n, a) + S cos 6* Jo' n (x/ jn , a) 



3tt 



(F[' n (xi tn , a) + S cos 9J l { n {xi tn , a)) 



(5.3) 



in the narrow width limit for the decay of the W boson. Expressions for Fq™ and Jg " can 
be found in [3]. Integrating over cos# gives us the differential energy distribution, 

dT i,n G 2^ m 5 - 



dxi, 



16vr 3 



Fo' n (xi n ,y) - -^-F l { n {xi n ,y) 



(5.4) 



We compared this theoretical prediction with the distribution obtained from the POWHEG 
method before interfacing with the Herwig++ parton shower. The best fit distributions 
shown in Figure 8 were obtained by setting as to 0.1 in (5.4). 




Figure 8: Scaled energy fractions of the charged lepton (left) and neutrino (right) from top decay. 
Black(solid)= Theory, Blue(dashes)= Decay. 



6. Truncated Shower 

The POWHEG method requires the addition of a 'truncated shower' before the hardest gluon 
emission in order to simulate the soft radiation distribution [21]. Due to angular ordering, 
the 'truncated' radiation is emitted at a wider angle than the angle of the hardest emission 
but at a lower pr- This means the 'truncated' radiation does not appreciably degrade 



- 16 - 



the energy entering the hardest emission and justifies our decision to generate the hardest 
emission first. 

In [18], there is a description of a method to generate a truncated shower of at most 
one gluon for the case of light quark production from e + e _ annihilation. In this section, 
we extend the discussion to top pair production. Below is an outline of how the 'truncated 
shower' was generated. We will consider the case in which at most one extra gluon is 
emitted by the top or anti-top before the hardest emission. The outline closely follows 
the Herwig++ parton shower evolution method described in [22, 23] where the evolution 
variables z, the momentum fractions, and q, the evolution scale, determine the kinematics 
of the shower. 

i) Having generated the pr of the hardest emission as discussed in Section 2 and the 
energy fractions x and y, calculate the light-cone momentum fractions z and 1 — z of 
the partons involved in the hardest emission. We will assume henceforth that x > y 
and that y is the energy fraction of the quark, i.e. the quark is involved in the hardest 
emission. Then 



Q-b + a g 



(6.1) 



where if we define 



m 



we have 



b = 

s 

A = \f{l — 46) (6.2) 

a(l + A) + y/x 2 (l + A) 2 - 8(6 + «)(! + A~^26) 
2(1 + A - 26) 



= y(l + A) - yg(i + A) 2 - 86(1 + A - 26) 
a ° 2(1 + A - 26) 

a g = - a h - a c (6.3) 

ii) Next generate the light-cone momentum fraction zt of the 'truncated' radiation within 
the range 

^ <Zt <l-9jL ( 6 .4) 

and distributed according to the massive splitting function, Pqq = Cf — Zt ^z t )q 2 

q~i is the initial evolution scale, i.e. y/s = 500 GeV, and Q g is a cutoff introduced to 
regularize soft gluon singularities in the splitting functions. In this report, a Q g value 
of 0.75 GeV was used, zt is the momentum fraction of the quark after emitting the 
'truncated' gluon with momentum fraction 1 — zt- 
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iii) Determine the scale of the hardest emission from 

/ Pt 2 m 2 



iv) Starting from an initial scale qi, the probability of there being an emission next at the 
scale q is given by 



S(qi,q)= (6.6) 



where 

A(q c ,q) = exp 



% / dz^p 00 e(o< 1 4.<p r ) • (6.7) 

q c is the lower cutoff of the parton shower which was set to the default value of 0.631 
GeV in this report, a s is the running coupling constant evaluated at z(l — z)q, Pqq 
is the Q — ► Qg splitting function and pt is the transverse momentum of the hardest 
emission. The Heaviside function ensures that the transverse momentum, p T of the 
truncated emission is real and is less than jyp. To evaluate the integral in (6.7), we 
overestimate the integrands and apply vetoes with weights as described in [22]. With 
r a random number between and 1, we then solve the equation 

S(q t ,q) = r (6.8) 

for q. If q > q^, the event has a 'truncated' emission. If q < q\ , there is no 'truncated' 
emission and the event is showered from the scale of the hardest emission. 



v) If there is a 'truncated' emission, the next step is to determine the transverse momen- 
tum p l T of the emission. This is given by 

p't = ^{l-z t f{z^-ml)-z t Q 2 . (6.9) 
If ptj? < or p f T > pt go to ii). 



vi) We now have values for zt, the momentum fraction of the quark after the first emission, 
p T , the transverse momentum of the first emission, z, the momentum fraction of the 
hardest emission and pr, the transverse momentum of the hardest emission. We can 
then reconstruct the momenta of the partons as described in [22]. The orientation of 
the quark, antiquark and hardest emission with respect to the beam axis is determined 
as explained there for the hard matrix element correction. 



In this paper, we consider only truncated emissions in the production process, not in the 
decay. 
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7. Parton shower distributions 



Next we interface the generated events with the Herwig++ 2.2.0 [24] parton shower and 
veto the hardest emissions in the production and decay of the top and anti-top pairs. In 
this section we will consider collisions at -^s = 500 GeV and only include the truncated 
shower for the production emissions. We considered four cases: 

1. Leading Order (LO): No POWHEG emissions. 

2. Production (Pr): Only POWHEG emissions in the production are allowed including the 
truncated shower. 

3. Decay (Dc): Only POWHEG emissions in the decays of the top/anti-top pairs are al- 
lowed. 

4. Production + Decay (PrDc): Both production and decay emissions are allowed. 

The following distributions were investigated in the lab frame for the two different e + e~ 
initial polarizations: 

i) The angle between the lepton from the decay of the top anti-quark and the top quark 
are presented in Figure 9. 

ii) The angle between the lepton and anti-lepton from the decays of the top pairs are 
presented in Figure 10. 

iii) The energy distributions of the b quark and b anti-quark before hadronization are 
presented in Figures 11 and 12. 

iv) The transverse momenta w.r.t the beam axis of the b quark and b anti-quark before 
hadronization are presented in Figures 13 and 14. 

v) The longitudinal momenta (along the beam axis) of the b quark and b anti-quark before 
hadronization are presented in Figures 15 and 16. 
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Figure 9: Angle between the lepton from the decay of the top anti-quark and the top quark. 



Angle between the charged lept oris from the top and anti-top decays Angle between the charged leptona from the top and anti-top decays 




8/pad 



Figure 10: Angle between the lepton and anti-lepton from the decays of the top pairs. 
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Figure 11: Energy of the b-quark before hadronization. 




Figure 12: Energy of the b anti-quark before hadronization. 



- 21 - 



Transverse momentum of b quark (GeV) Transverse momentum of b quark (GeV) 




Figure 13: Transverse momentum of the b quark before hadronization. 



Transverse momentum of b anti— quark (GeV) Transverse momentum of b anti— quark (GeV) 




Figure 14: Transverse momentum of the b anti-quark before hadronization. 
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Figure 15: Longitudinal momentum of the b quark before hadronization. 




Figure 16: Longitudinal momentum of the b anti-quark before hadronization. 
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At leading order, the leptonic correlations in Figures 9 and 10 are as expected with 
higher correlations seen for e^e^ annihilation than for e^e^ annihilation. At next-to- 
leading order, it can be observed that the POWHEG emissions do not change the shapes of 
the distributions much except for a slight broadening of the peaks. 

Also at leading order, the distributions in Figures 11-16 have the expected shapes with 
the b quarks and anti-quarks having softer (harder) energy, longitudinal momentum and 
transverse momentum spectra for e^e^ (e^e^,) annihilation. 

Now comparing the POWHEG production and decay distributions for the b quark in 
Figures 11-16, we observe that the decay emissions soften the spectra more than the pro- 
duction emissions and therefore these have the greater effect in the production + decay 
distributions. This is expected since the scale range available for the production emissions 
« \og{-\fs/m\) is less than the range available for the decay emissions « log(m t /mb). 

8. Conclusions 

Using the Monte Carlo event generator Herwig++, we have successfully applied the POWHEG 
method to investigate angular correlation distributions at next-to-leading order in top pair 
production and decays at ILC energies. In all distributions studied, the POWHEG emissions 
have the effect of broadening the peaks of the leading order predictions slightly. We also 
compared momentum distributions of the b quarks and anti-quarks before hadronization 
and observe that the decay emissions soften the spectra more at next-to-leading order as 
expected. 
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Abstract: We study the effects of gluon radiation in top pair production and their decays 
for e + e~ annihilation at the ILC. To achieve this we apply the POWHEG method and interface 
our results to the Monte Carlo event generator Herwig++. We consider a center of mass 
energy of y/s = 500 GeV and compare decay correlations and bottom quark and anti-quark 
distributions before hadronization. 
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1. Introduction 

High energy polarized e + e~ colliders will be essential instruments in the search for the fun- 
damental constituents of matter and their interactions. One such collider being designed is 
the International Linear Collider (ILC) which is expected to run at centre-of-mass energies 
> 500 GeV. At these energies there will be a significant proportion of top quark pairs 
produced from the annihilation process so the ILC provides an impressive tool to carry out 
detailed studies of top quark physics. The top quark radiates gluons both in its produc- 
tion phase and its decay phase, thus it is useful to see how the leading order experimental 
analysis will be affected by these QCD corrections. There have been numerous studies of 
top quark production and their decays both at leading order and next-to leading order in 
QCD. Details can be found in [1-8] and many more besides. In this paper we consider 
the process at next-to-leading order in the production and semi-leptonic decays of the top 
pairs using the POWHEG method [9, 10] in conjunction with the Monte Carlo event generator 
Herwig++ [11]. The POWHEG method has been successfully applied to Z pair production [12], 
heavy flavour production [13], e + e~ annihilation into hadrons and Drell-Yan vector boson 
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production [14, 15]. We work in the narrow width approximation and hence do not in- 
clude interference between the production and decay emissions which are negligible in this 
limit [16, 17]. We also take account of the beam polarization and spin correlations of the 
top pairs. Finally, we present plots of some relevant distributions. 



2. Hardest emission generation: Production 



The order-a*. differential cross section for the process e + e 
a vector current, can be written as 



V — > ttg where V represents 



R{x,y) = <7 V W v (x,y) 



oy 2a s 

V 37T 



(x + 2pf + (y + 2pf + Cv 



2p 



2p 



(l + 2p)(l-x)(l-y) (l- x ) 2 (l-yf 

(2.1) 

where ay is the Born cross section for heavy quark production, Wy = Rjoy, p = rn^/s 
where m t is the mass of the top quark and s is the square of the center of mass energy, 
Qy = — 8/?(l + 2p), v = ^Jl — 4p and x, y are the energy fractions of t and i respectively. 
In the case of the axial current contribution e + e~ — > A — > ttg, we have 



R{x,y) = a A W A {x,y) 



v 3ir 



(x + 2pf + (y + 2pf + Ca 



2p 



2p 



(1 - Ap)(l - x)(l - y) (l- x f (\- y y 



(2.2) 

where a a is the Born cross section for heavy quark production by the axial current, Wa = 
RjoA and (a = 2p[(3 + x g ) 2 — 19 + 4p] where x g = 2 — x — y. 

Because of the top mass, the phase space for gluon emission is reduced and the collinear 
divergences present in the massless quark cross-section are regularized here. However, the 
infra-red divergence as the gluon momentum goes to zero is still present. We can write 
down the cross-section for the hardest emission as 



do = ^2,B{v)d® % 



1 R 



B(v) 



(2.3) 



where B{v) is the Born cross section and v represents the Born variables, r represents 
the radiation variables and dQ v and d$ r are the Born and real emission phase spaces 
respectively. A^ lo (pt) is the modified Sudakov form factor for the hardest emission with 
transverse momentum px, as indicated by the Heaviside function in the exponent of (2.4), 



A™( PT ) = exp 



/ d$r^^-e(k T (v,r)-p T 



Furthermore, 



B(v) = B{v) + V(v) + j {R{v, r) - C{v, r))d<5> r . 



(2.4) 



(2.5) 



B(v) is the sum of the Born, B(v), virtual, V{v) and real, R(v,r) terms, (with some 
counter-terms, C(v,r)). It overcomes the problem of negative weights since in the region 
where B(v) is negative, the NLO negative terms must have overcome the Born term and 
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hence perturbation theory must have failed. It is used to generate the variables of the Born 
subprocess to which the real-emission contributions factorize in the collinear limit. 
Now explicitly for e + e~ — > ttg, 



A^°(p T ) = exp 



J dxdyW{x,y)Q(k T (x,y)) - pr) 



where 



k T (x,y) 



[1 - x)(l -y)(x + y-l)-p(2-x- y) 2 
max(a;, y) 2 — 4p 



(2.6) 



(2.7) 



is the transverse momentum of the hardest emitted gluon relative to the splitting axis, as 
illustrated in Figure 1 below. 




Figure 1: Transverse momentum, kx- 



2.1 Generation of radiation variables, x and y 

The radiation variables, x and y are to be generated according to the probability distribu- 
tion 

& w {k T )W{x,y)dxdy (2.8) 

where for e + e~ annihilation via a vector and axial current, W(x, y) and A w (kr) are defined 
in (2.1), (2.2) and (2.6). 

In the region where x > y, let us define the dimensionless variable, k as 



K = 4 = (I ~ - y)(x + y - 1) - p(2 - x - yf 
s x 2 — 4p 

There are two solutions for y for each value of x and k. 



(2.9) 



2/1,2 



3a: - 2px + 2 + 4p ± ^{x 2 - Ap){An{x - 1 - p) + (x - l) 2 ) 
2(1 + p-x) 



(2.10) 



3 



/ 



Exchanging the y variable for k, we find 

W{x,y)@{k T {x,y) - p T )dxdy = I dx I dn — s - —W(x,yi 2 ) 

Jx mi „ Jk 3tt dn 



I dX I 



2a s (K.s) \fx 2 - 4p 

dn o r, A , = ; , 1/1,2) 

3tt ^/(4K(a; - 1 - p) - (a; - l) 2 ) 

(2.11) 



Now the integrand, W' = V x 4 p =W(a:,j/i2) in (2.11) yields a complicated 

y(4/t(x-l-p)-(x-l) 2 ) 

integral so we look for an upper bound on w' which we denote as V' = ^V{%, 2/1,2) to 
simplify the integration. We then generate the radiation variables as outlined below: 

1. Set p m ax = ^Tmax- 



2. For a random number, n between and 1, solve the equation below for pt 

A V ( PT ) 

n = A^y (2 - 12) 

where A v (pt) = exp — J dx dnV' (x, k) 

3. Generate the variables x and y according to the distribution 

V{x,y)6{hr{x,y)-pT). (2.13) 

4. Accept the generated value of pr with probability W/V. If the event is rejected set 
Pmax =PT and go to 2). 

Using our knowledge of the form of the integrand in the massless quark case [18], we 
guess that V' should take the form, 

V(x,k)=N k sK — — 2.14 

37r (1 — x + 7(k, a;))7(K, x) 

where 

7(k, a:) = y/ (1 - x)(l - 2k - 2^?TpK) . (2.15) 

and A^ K is a normalisation factor which depends on k, which has to be tuned to ensure 
that V' is an upper bound of W' . Both V' and W' have the same divergent behaviour at 

a^max = 1 - 2K - 2a/« 2 + pK. 

Table 1 shows the N K values for different ranges of k and the two solutions for y used for 
both axial and vector currents. The lower limit on k was set by choosing kx = Aqcd = 0.2 
GeV thus setting a lower bound on the transverse momentum. 

We now consider the specific case where ra t = 175 GeV and y/s = 500 GeV i.e. 
p = 0.1225. For k < 0.024, there are two y solutions in the region of phase space where 
x > y. This is illustrated in Figure 2 below for k = 0.01. The red line denotes the phase 
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Range of k 




N K (y 2 ) 


0.024 - 0.03 


0.4 


0.4 


0.015 - 0.024 


0.7 


0.7 


0.005 - 0.015 


1.2 


1.1 


0.0005 - 0.005 


4.0 


2.6 


0.0001 - 0.0005 


9.0 


6.0 


0.00005 - 0.0001 


13.0 


7.0 


0.00003 - 0.00005 


17.0 


10.0 


0.00000016 - 0.00003 


45.0 


35.0 



Table 1 : N K for different values of n for both axial and vector currents 




Figure 2: Phase space and y solutions for k = 0.01 in the region x > y. 



space for gluon emissions. The two solutions lie on either side of the dashed line 

(2-x)(l-x) + 2p) 



2(1 -x + p) 



(2.16) 



and are equal when 



1 - 2k - 2a/k 2 + pn (2.17) 



which lies on the dotted line. At n = 0.024, the branches meet along the line y = x and 
there is only one solution for y in the region (the lower branch). So for k > 0.024, only 
one y solution exists for x > y. This is illustrated in Figure 3 for n = 0.028 . In addition 
there are no y solutions for k > 0.03 for x > y . Also note that for x < x u , there is only 
one solution for y. 
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Figure 3: Phase space and y solutions for k = 0.028 in the region x > y. 

In the region where there are two solutions, the integral (with V in place of W ) in 
(2.11) is performed along both branches independently and summed. For the upper branch, 



x runs from x u to £ max = 1 — 2k — 2\J + pn while for the lower branch, x runs from x 
to s max where if we define 



x a = 39k - 1 + 12p - 168pn - 48p z + k a + lhn 1 + \2pn l + 48/cTk + Up 6 , 
x b = 6(-33k 2 + 3k - 3k 3 + 288p 2 K - 768/9 3 k - 48pn + 168,0k 3 - 204p 2 K 2 + 300pK 2 
+ 12pK 4 + 768p 4 K + 144p 2 K 3 + 486p 3 K 2 )^ , 



X c — "s/^fl 2 %b 
, -1 ( x b 

Xd = tan — 



1 I f x A 

Xr COS 

12 V 3 J 

(-1 - k 2 - 10k + 8p - 2pn - 16p 2 ) cos 

12x1 

K + 5 + 4p 
6 ' 

V^3 fx d \ f j , 1 + k 2 + 10k - 8p + 2pK + 16p 2 

"X" 

,3 



Xh 



v| sin ^ r§ + i+- 2 +^-y+2 P K + i6p 2 5 (2 18) 



we can write a;" and x as 



3v — X g | 3/ f | ^ | 3/ fa . 



£ e + Xf + a; 9 - a/, ; ( 2 -19) 
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In the region where there is only one solution for y, x runs from x l to x u . 
The k integration can then be performed numerically. Having performed the integration, 
values for k and hence kx are then generated according to steps 1 and 2 in Section 2.1. 
The variables x and y are then to be distributed according to W'(x, y)5(kT{x, y) — pr)- 
This is the subject of the next section. 

2.2 Distributing x and y according to W(x,y) 

To generate x and y values with a distribution proportional to V(x,y)S(kT(x,y) — pr), 
where from (2.11), V(x,y) is the d ^ dK , we can use the 8 function to eliminate the y 
variable by computing 

V{x,y) 



D(x) = I dyS(k T -p T )V{x,y) 



dk-. 



(2.20) 



where y is such that kr{x,y) = pr- Note that ^ 1 is the same for both y solutions. We 
then generate x values with a probability distribution proportional to D with hit-and-miss 
techniques as described below. All events generated have uniform weights. 

1. Randomly sample x, N x times (we used N x = 10 5 ) in the range [x m i Q : x max ] for the 
selected value of k. 



2. For each value of x, evaluate D = D(x, y\) + D(x, y 2 ) if there are 2 solutions for the 
selected k and D = D(x,y2) if there is only one solution. Also, if k < 0.024 and 
x < x u (see Figure 2), there is only one y solution so evaluate D = D(x,y2)- 

3. Find the maximum value D max of D for the selected value of k. 

4. Next, select a value for x in the allowed range and evaluate D. 

5. If D > rD max (where r is a random number between and 1), accept the event, 
otherwise go to 4.) and generate a new value for x. 

6. If for the chosen value of x, there are two solutions for y, select a value for y in the 
ratio D(x,yi) : D(x,y 2 ). 

7. Compare V' (x,y) with the true integrand, W'(x,y). If the event fails this veto, set 
K max = k and regenerate a new k value as discussed in Section 2.1. 



NB: For the region y > x, exchange x and y in the above discussion. In this way, the 
smooth phase space distribution in Figure 5 below was obtained for the hardest emission 
events for an axial current. The plot show 2,500 of these events. 
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Figure 4: Phase space and distribution of hardest emissions for axial (left) and vector(right) 
currents with p = 0.1225. 
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Figure 5: Phase space and distribution of hardest emissions for an axial current with p = 0.0625 
(left) and p = 0.01361 (right). 



The procedure was repeated for y/s = 700 GeV, p = 0.0625 and y/s = 1500 GeV, 
p = 0.01361 and the corresponding plots are shown below. As can be seen, the method is 
stable as p — > 0. This is not surprising because the upper bound function V(x, k) in (2.14) 
is stable as p — > and tends to W' , the true value of the integrand in this limit. 
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3. Hardest emission generation: Decays 



In addition to gluon emission in top production, we also studied the emission in its decay, 

t(pi) -> W + (wi)b(ri)g(k) . (3.1) 

The procedure for generating the hardest emission in this case follows the same lines as 
discussed in Section 2. We parameterize the phase space for the decay in terms of variables 
x and y defined as 



2w\ ■ pi 



y - 2 



x = ^ (3.2) 

m 2 

where a = m^/m 2 with m w and m t the masses of the W boson and top quark respectively. 
(y + a)/2 and x/2 are the energy fractions of the W boson and gluon in the top frame. 
Therefore the corresponding energy fraction of the b quark in this frame is given by 

Xh 2 — y — a — x .„ „. 

f = —2 • (3 - 3) 

In this paper, we neglect the b mass and work in the narrow-width approximation so that 
the top quarks and W boson are on-shell. The t — > Wbg differential decay rate is given by: 



1 d 2 T a s C F 



Tq dxdy tt (1 — y)x 2 



(1 - y){\ - x) + x 2 (y + x-1) 2 2a{\-y)x 



+ X—— tt, h 



2(1 -a) 2 (1 -a) 2 (l + 2a) 

(3.4) 

where Tq is the leading order decay rate. The phase space limits for the decay are: 

+ {l-x)<y<l, 



i - x 

0<x<l-a. (3.5) 

Working in the rest frame of the top quark where the parton shower is formulated in 
Herwig++, we identify the splitting axis corresponding to the original b — W boson axis 
and therefore the relative transverse momentum for gluon emission is: 



, , x / (I ~y)(y + x(2-y-a)- x 2 -l) 

k T (x,y)=m i] j . (3.6) 

k 2 

Now defining a dimensionless variable k = -jfe, we find that in analogy to the production 
case, there are 2 solutions for y for each value of x and k. 

_ x 2 +ax + 2-3x- 2an ± ^J(x 2 - 4k(1 + a))(x - l) 2 + Aan(An + 1 - a) + x 2 (a + 2x - 2) 
yi ' 2 = 2{k + 1-x) 

(3.7) 

These solutions may be identified with either initial state gluon emission from the top 
quark (j/2) or final state radiation from the bottom quark (y\). A plot of the phase space 
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Figure 6: Phase space(solid), y (dot-dash) and solutions yi (dots) and yi (dashes). 



and the 2 solutions for k = 0.01 is shown in Figure 6. We then construct the modified 
Sudakov form factor for the generation of the hardest emission. The exponent of the form 
factor is given by 

[ W(x : y)@(k T (x,y) -p T )dxdy = f^" dx f^" dK ^™^L ^-W(x, k) , (3.8) 
J ix min Jk n dK 

where W(x,n) is the differential cross-section (3.4) and ^| is the Jacobian for the change 
of variables from y to k. Note that k = K max = ^— when the W boson is at rest and 
x = 1 — ^/a, y = For a given k, a; m ; n and x max are also given by 

Xm&x = 1 - a - . (3.9) 

To make the integral simpler, we again look for an upper bound V (x, k) on the integrand 
as we did for the production case. To do this we replace the Jacobian with the simpler 
expression, 

dy d f x 2 — 3x — 2na + 2 + xa\ —a x 2 — 3x — 2na + 2 + xa 

dK~dn\ 2(k + 1-x) )~k+1-x 2(k + 1 - x) 2 

where y lies in between y\ and yi and is indicated in Figure 6. We also overestimate the 
differential cross-section by replacing (3.4) with 

asCp I — a 



(3.10) 



U(x,y) =N K - 



7T 2z 2 (l-y') 



(3.11) 



where N K is a normalisation factor dependent on k and is chosen such that V' = U^j- is 



greater than the integrand in (3.8). The N K values are given in Table 2 for the 2 solutions. 
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Range of k 



iV K ( yi )xl0 4 iV K (y 2 )xl0 4 



a a i a A'70'7 

0.01 — 0.0737 


0.005 


0.006 


0.005 — 0.01 


0.0175 


0.02 


A f~\f~\1 A A A P" 

0.001 — 0.005 


A AO 

0.03 


A A 1 P" 

0.045 


A AAA1 A AA1 

0.0001 — 0.001 


A AO 

0.08 


A 1 A 

0.12 


A AAAAT A AAA1 

0.00005 — 0.0001 


A A 

0.2 


A A 

0.2 


A AAAAO f A AAAAr" 

0.000025 — 0.00005 


A O 

0.3 


A A 

0.2 


0.0000075 - 0.000025 


1.0 


0.9 


0.000005 - 0.0000075 


2.0 


0.9 


0.0000025 - 0.000005 


3.0 


0.9 


0.0000013 - 0.0000025 


6.0 


0.9 



Table 2: N K for different values of k 

1 



0.95 
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Figure 7: Phase space distribution of P0WHEG events 

The lower limit on k is 1.3 x 10~ 6 and was set by choosing kr = Aqcd = 0.2 GeV thus 
setting a lower bound on the transverse momentum. We then generate the values of k and 
distribute x and y according to the true differential (3.4) using vetoes as described for the 
production case in Sections 2. Figure 7 shows the phase space distribution obtained. 

4. Spin Correlations and the distribution of Born variables 

In [19], it was observed that the lepton matrix element for the production process 

e + (p) + e-(q) -> t( Pl ) + t(p 2 ) + g(p 3 ) -> W + {wi) + 6(n) + W-{w 2 ) + b(r 2 ) + g(p 3 ) 

-> l + {k-i) + u{x x ) + b{n) +r(k 2 ) + v{x 2 ) + b(r 2 ) + g(ps) (4.1) 
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is bounded from above in the narrow width approximation by the undecayed matrix ele- 
ment obtained by eliminating the decay products i.e. W + , W~ , 6, b and putting the parent 
particles i.e. t, i on-shell, multiplied by a process dependent constant. We can then use the 
undecayed matrix elements to perform computer-intensive tasks such as event generation 
and finally, by using the hit-and-miss method, replace the parent particles with their decay 
products. This procedure is outlined below: 

1. Evaluate the undecayed matrix elements which are proportional to the upper bounds 
on the lepton matrix elements. Generate hard events using the POWHEG method 
described above with the top and anti-top quarks in the final state. 

2. For each event, generate the decay products and their four-momenta according to 
the phase space. 

3. Evaluate the leptonic decay matrix element for each event. If the decay matrix 
element divided by the corresponding upper bound is less than a random number r 
between and 1, throw away the decay momenta and return to step 2. 

4. Otherwise, replace the top and anti-top momenta with the decay momenta and shower 
the event. 



4.1 Undecayed matrix elements 

At the ILC, the electron and positron beams will be polarized i.e. either e^e^ or e^e^ 
where the subscripts L and R represent the left-handed and right-handed helicity states 
respectively. The corresponding undecayed matrix element for e^e^ annihilation is: 



M(e L (p)e+{q) -> t St { Pl )i Sl {p 2 )g{p 3 )) 
[v{q)7»Lu{p)]u{pi,s t ) — i ( a " 



_2p 2 -p 3 V 2 

1 / - . - . \ / a LL n . O-LR n 

+ ^v~3 lv{vi +n+ mt) y~ lL + ~ lR 



7l + ^tIr) {-P2 ~ P3 + mthv 



T a v(p2,StK(p 3 ) 



(4.2) 



where e is the polarization vector of the gluon, T a is the colour matrix, p = p^7^ and 
7r/l = 7^(1 i 75) /2- For e^e\ annihilation, interchange L, R in the above equation, st, 
are the spin vectors of the top and anti-top quarks respectively and satisfy the relations: 



s t ■ p\ = 

St ■ P2 = 
Sf s t = -1 
Si ■ st = -1 



(4.3) 
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The massive spinors u(p, s), v(p, s) are given in terms of the massless spinors u(p),v(p) by 



«(P,t) = 2 

u (p,i) = 2 75S ^(p) 
«(P,t) = 3 «(?) 

«(P)4.) = n 75 ^ (p) 



(4.4) 



The coupling constants au are given by 

Qt + QlQt "V 2 



e 2 g 

au = 



sin 2 W s - Ml + «M z r z 



(4.5) 



where Mz is the Z boson mass, Y z is the width of the Z boson, Qyy is the Weinberg angle, 
Qt is the electric charge of the top in units of the electric charge e, g = ^4nas and s is 
the center of mass energy squared. The couplings to the Z boson are given by 

L _ 2 sin 2 6 W - 1 



Q 



2 cos Qw 
R sin 2 6>vf 



Q 



cos 

L _ 3 - 4 sin 2 6 W 



6 cos #vk 



= 

3 cos 6% 

In Section 2.2, we distributed our events according to the vector and axial vector 
current matrix elements separately using the POWHEG method. To obtain a full unpolar- 
ized distribution we can select events from either current distribution according to their 
contributions to the full cross-section given below. 

a = 3/3(1 + 2p) (l + ci^) a vv + 3/3 3 (l + di^) o AA 

avv = [ Q 2 _ 2 Q t V e V tX i(s) + (A 2 + V e 2 )V t 2 X2 (s)} 

s 

-AA = ^[(A 2 e + V 2 )A 2 tX2 (s)] , (4.7) 
s 

where (3 = ^/l — p, a em is the electromagnetic coupling, At, A e and Vt, V e are the axial and 
vector coupling constants of the top t and electron e to the Z boson and c\ = 3.5 and 
d\ = 2.25 are the QCD correction coefficients defined at m t = 175 GeV and ^/s = 500 GeV 
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i.e. p = 0.1225 [20]. Xi( s ) an d X2{s) are given by 



Xi(s) 



'{s-MlY + YlMl 

X2 ( S ) = = K 



2 £_ 



V2G F M 2 7 



167ra em 



(4.8) 



where Gf is the Fermi constant and Mz and Y z are the mass and decay width of the Z 
boson respectively. 

Explicit expressions for the Born, virtual and real polarization dependent squared 
matrix elements for the production process are given in [5] in terms of the energy fractions 
x, y of the top and anti-top quarks and the polar angle and azimuthal angles orienting 
the tig plane relative to the e + e~ beam axis. For each initial polarization, we then assign 
final-state polarizations to each event in proportion to the squared matrix elements and 
distribute the polar and azimuthal angles of the top/anti-top pairs accordingly using well- 
known Monte Carlo techniques. 

4.2 Decay matrix elements 

Next, we investigate the decays of the top and anti-top pair. The leptonic matrix elements 
for the process in (4.1) are dependent on the spins of the top and anti-top quark. This 
dependence can be written in the form of a decay density matrix. The decay density matrix 
p x A ' , for an on-shell top quark is given by 

p ^ 

K - mlY + {m w T w y 

(r-i • xi){pi ■ ki) - (s L ■ &i)(ri • xi)m t -{k x ■ n)(x 1 ■ ri)m t - ie(pi, hi, s t , n)(x 1 ■ r x ) 
-(ki ■ n)(xi ■ ri)m t + ie(p 1 ,k 1 ,s t ,n)(x 1 ■ n) (n • xi)(pi ■ hi) + (s t ■ ki)(n ■ xi)m t 

where A, A' are spin labels, st is the top spin vector, n is a spacelike vector perpendicular to 
st and pi and m w , Ty/ are the mass and width of the W boson respectively. In this paper 
we work in the helicity basis for which the top quark spin is defined along its direction 
of motion. A similar matrix can be derived for i decay. The spin-specific decay matrix 
elements are therefore of the form: 

= (4-9) 

where At, Aj are spin labels for the top and anti-top respectively and M is the matrix 
element for the undecayed process introduced in Section 4.1. By diagonalizing the density 
matrix, we can obtain the largest possible value of the matrix elements and hence the upper 
bound. An explicit computation gives this upper bound | M* 6 | 2 on the top decay as [19], 

|M* l 2 = | Vw I 2 (ri_-£iKgi_ L fci) . ^ |2 . . 

1 ub 1 [K - ml )2 + {m w T w f}[{p\ - m\f - (m t T t f] 11 1 ' ' 
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where | M | 2 is the undecayed matrix element for unpolarized tig production. A similar 
expression | M* 6 | 2 can be obtained for the decay of the top anti-quark by interchanging 
the labels 1 and 2 and t and i in (4.10). Hence the full upper bound can be written as: 

I M t l 2 l M* I 2 

| M tt |2 = I M ub I J M ub I (4n) 

| M | 2 

Having obtained the decay matrix elements and their upper bounds, we then proceed 
to generate events with leptons in the final state as outlined at the beginning of this section. 

For the POWHEG decays, we apply the same method where in this case the undecayed 
matrix elements \M\ are the leading order matrix elements for the process 

e+ + e~ -+t + i . (4.12) 



We then use the next-to-leading order decay matrix for which the helicity amplitudes can 
be found in [5]. These are given in terms of the polar and azimuthal angles of the decay 
w.r.t the top/anti-top axis and we distribute them as described for the production process 
in Section 4.1. Note that in this case, we generate two decay gluons, one each from the top 
and anti-top quark. 

In addition, we also consider POWHEG radiations in both the production and decay 
process by independently generating the emission and distributing the Born variables of 
the production process first and then generating the emission and distributing the Born 
variables of the decay process to yield three gluons in the final state. 



5. Decay NLO lepton spectra comparisons 

Extensive studies have been carried out on the lepton angular and energy distributions 
from the semi-leptonic decays of polarized top and anti-top quarks at next-to-leading order 
in a s [1,3]. 

t W + + b + g^re + + v e + b + g 
i -+W~ + b + g -s- e" + u e + b + g . (5.1) 

In the top rest frame, we define 6 as the angle between the spin 3- vectors St,sj of the 
decaying quark and the lepton. We also defined the scaled energies x l,n of the charged 
lepton and the neutrino respectively as 

2E t 

xi = 

m t 

x n = (5.2) 

m t 

where E\ and E n are the energies of the charged lepton and neutrino in the top rest 
frame. In these variables, the NLO double differential distribution of the charged lepton 
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and neutrino in the decay of a heavy top or anti-top quark with polarization S has been 
shown to be of the form 

dT l > n _ G 2 F ml 



dx l > n d cos 32vr 3 



Fi> n (x l > n ,a) + Scos6jk n (x l * n ,a) 



^(Fj'V' n ,«) + Scosfljj'V'V)) 

07T 



(5.3) 



in the narrow width limit for the decay of the W boson. Expressions for Fq* and Jq" can 
be found in [3]. Integrating over cos 6 gives us the differential energy distribution, 



Ln 



dx l ' n 



16tt 3 



o 



x l > n ,y) 



2as r>l,n,l 



3tt 



(5.4) 



We compared this theoretical prediction with the distribution obtained from the POWHEG 
method before interfacing with the Herwig++ parton shower. The best fit distributions 
shown in Figure 8 were obtained by setting as to 0.1 in (5.4). 




Figure 8: Scaled energy fractions of the charged lepton (left) and neutrino (right) from top decay. 
Black(solid)= Theory, Blue(dashes)= Decay. 



6. Truncated Shower 

The POWHEG method requires the addition of a 'truncated shower' before the hardest gluon 
emission in order to simulate the soft radiation distribution [21]. Due to angular ordering, 
the 'truncated' radiation is emitted at a wider angle than the angle of the hardest emission 
but at a lower py. This means the 'truncated' radiation does not appreciably degrade 
the energy entering the hardest emission and justifies our decision to generate the hardest 
emission first. 
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In [18], there is a description of a method to generate a truncated shower of at most 
one gluon for the case of light quark production from e + e~ annihilation. In this section, 
we extend the discussion to top pair production. Below is an outline of how the 'truncated 
shower' was generated. We will consider the case in which at most one extra gluon is 
emitted by the top or anti-top before the hardest emission. The outline closely follows 
the Herwig++ parton shower evolution method described in [22, 23] where the evolution 
variables z, the momentum fractions, and q, the evolution scale, determine the kinematics 
of the shower. 

i) Having generated the pr of the hardest emission as discussed in Section 2 and the 
energy fractions x and y, calculate the light-cone momentum fractions z and 1 — z of 
the partons involved in the hardest emission. We will assume henceforth that x > y 
and that y is the energy fraction of the quark, i.e. the quark is involved in the hardest 
emission. Then 

(6.1) 



otb + en, 



where if we define 



6 = ^ 



A = \f{\ - 46) (6.2) 



we have 



otr. 



x(l + A) + y/x 2 {l + A) 2 - 8(6 + k)(1 + A - 2b) 
2(1 + A - 26) 

y(l + A) - y/g(I + A) 2 - 86(1 + A - 26) 
2(1 + A- 26) 

2 

j- — ^ - a b -a c (6.3) 

(6.4) 



ii) Next generate the light-cone momentum fraction zt of the 'truncated' radiation within 
the range 

— < z t < 1 — -T- (6.5) 



l+z t 2 2m\ 



l-z t z t (l-z t )q 2 



and distributed according to the massive splitting function, Pqq = Cf 
qi is the initial evolution scale, i.e. y/s = 500 GeV, and Q g is a cutoff introduced to 
regularize soft gluon singularities in the splitting functions. In this report, a Q g value 
of 0.75 GeV was used, zt is the momentum fraction of the quark after emitting the 
'truncated' gluon with momentum fraction 1 — zt- 



iii) Determine the scale q^ of the hardest emission from 

_ / PT 2 ~ ml Q g 2 

Qh ~ \ ~2M ^ + T2" + ^2 

y z z (l — z) z z(l — z) 
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iv) Starting from an initial scale the probability of there being an emission next at the 
scale q is given by 

A(<? c ,<?j) 



where 



A(<Zc,<z) = exp 



S{qi,q) 



he ? J 



dz^-p QQ e(o < P f T < pt) 



(6.7) 



(6.8) 



q c is the lower cutoff of the parton shower which was set to the default value of 0.631 
GeV in this report, a s is the running coupling constant evaluated at z(l — z)q, Pqq 
is the Q — > Qg splitting function and is the transverse momentum of the hardest 
emission. The Heaviside function ensures that the transverse momentum, pip of the 
truncated emission is real and is less than px- To evaluate the integral in (6.8), we 
overestimate the integrands and apply vetoes with weights as described in [22]. With 
r a random number between and 1, we then solve the equation 



S(qi,q) = r 



(6.9) 



for q. If q > q^ the event has a 'truncated' emission. If q < q^ , there is no 'truncated' 
emission and the event is showered from the scale of the hardest emission. 



v) If there is a 'truncated' emission, the next step is to determine the transverse momen- 
tum pip of the emission. This is given by 

pip = ^(l-z t )Hz?<f - m \)-z t Q? . (6.10) 

If pip 2 < or pip > pt go to ii) . 



vi) We now have values for zt, the momentum fraction of the quark after the first emission, 
pip, the transverse momentum of the first emission, z, the momentum fraction of the 
hardest emission and pr, the transverse momentum of the hardest emission. We can 
then reconstruct the momenta of the partons as described in [22]. The orientation of 
the quark, antiquark and hardest emission with respect to the beam axis is determined 
as explained there for the hard matrix element correction. 



In this paper, we consider only truncated emissions in the production process, not in the 
decay. 
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7. Parton shower distributions 



Next we interface the generated events with the Herwig++ 2.2.0 [24] parton shower and 
veto the hardest emissions in the production and decay of the top and anti-top pairs. In 
this section we will consider collisions at ^fs = 500 GeV and only include the truncated 
shower for the production emissions. We considered four cases: 

1. Leading Order (LO): No POWHEG emissions. 

2. Production (Pr): Only POWHEG emissions in the production are allowed including the 
truncated shower. 

3. Decay (Dc): Only POWHEG emissions in the decays of the top/anti-top pairs are al- 
lowed. 

4. Production + Decay (PrDc): Both production and decay emissions are allowed. 

The following distributions were investigated in the lab frame for the two different e + e~ 
initial polarizations: 

i) The angle between the lepton from the decay of the top anti-quark and the top quark 
are presented in Figure 9. 

ii) The angle between the lepton and anti-lepton from the decays of the top pairs are 
presented in Figure 10. 

iii) The energy distributions of the b quark and b anti-quark before hadronization are 
presented in Figures 11 and 12. 

iv) The transverse momenta w.r.t the beam axis of the b quark and b anti-quark before 
hadronization are presented in Figures 13 and 14. 

v) The longitudinal momenta (along the beam axis) of the b quark and b anti-quark before 
hadronization are presented in Figures 15 and 16. 
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Figure 9: Angle between the lepton from the decay of the top anti-quark and the top quark. 




Figure 10: Angle between the lepton and anti-lepton from the decays of the top pairs. 
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Figure 11: Energy of the b-quark before hadronization. 




Figure 12: Energy of the b anti-quark before hadronization. 
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Transverse momentum of b quark (GeV) Transverse momentum of b quark (GeV) 




Figure 13: Transverse momentum of the b quark before hadronization. 



Transverse momentum of b anti — quark (GeV) Transverse momentum of b anti— quark (GeV) 




Figure 14: Transverse momentum of the b anti-quark before hadronization. 
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Figure 15: Longitudinal momentum of the b quark before hadronization. 




Figure 16: Longitudinal momentum of the b anti-quark before hadronization. 
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At leading order, the leptonic correlations in Figures 9 and 10 are as expected with 
higher correlations seen for e^e^ annihilation than for e^e^ annihilation. At next-to- 
leading order, it can be observed that the POWHEG emissions do not change the shapes of 
the distributions much except for a slight broadening of the peaks. 

Also at leading order, the distributions in Figures 11-16 have the expected shapes with 
the b quarks and anti-quarks having softer (harder) energy, longitudinal momentum and 
transverse momentum spectra for e^e^ ( e £ e #) annihilation. 

Now comparing the POWHEG production and decay distributions for the b quark in 
Figures 11-16, we observe that the decay emissions soften the spectra more than the pro- 
duction emissions and therefore these have the greater effect in the production + decay 
distributions. This is expected since the scale range available for the production emissions 
~ \og{sfs I mt) is less than the range available for the decay emissions ps log(m t /mb). 

8. Conclusions 

Using the Monte Carlo event generator Herwig++, we have successfully applied the POWHEG 
method to investigate angular correlation distributions at next-to-leading order in top pair 
production and decays at ILC energies. In all distributions studied, the POWHEG emissions 
have the effect of broadening the peaks of the leading order predictions slightly. We also 
compared momentum distributions of the b quarks and anti-quarks before hadronization 
and observe that the decay emissions soften the spectra more at next-to-leading order as 
expected. 
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